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Statistical  normalization  of 
non-Rayleigh  reverberation 

D.  A.  Abraham 


Executive  Summary:  The  primary  difficulty  in  the  detection  of  sub¬ 

marines  using  low-frequency  active  sonar  in  shallow  water  is  reverberation  and 
the  numerous  false  alarms  it  produces.  Automatic  detection  algorithms  are 
required  to  analyze  this  plethora  of  alarms  as  well  as  to  reduce  operator  load¬ 
ing  and  to  improve  consistency  in  performance.  Hindering  the  implementation 
of  automatic  detection  algorithms  is  the  variability  of  the  statistics  of  the  re¬ 
verberation  from  the  traditionally  assumed  Rayleigh  distribution,  resulting  in 
increased  and  unknown  false  alarm  performance  that  changes  in  range  and 
bearing.  This  report  proffers  an  algorithm  capable  of  removing  the  statisti¬ 
cal  variation  of  the  background  reverberation,  allowing  the  implementation  of 
automatic  detection  algorithms. 
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Statistical  normalization  of 
non-Rayleigh  reverberation 

D.  A.  Abraham 


Abstract:  Low-frequency  active  sonar  systems  operating  in  shallow  wa¬ 

ter  are  primarily  limited  by  reverberation.  Reverberation  is  traditionally  as¬ 
sumed  to  follow  a  Rayleigh  probability  distribution,  from  which  the  detector 
and  other  signal  processing  algorithms  are  developed.  Experimental  studies 
have  shown  that  reverberation  can  be  non-Rayleigh  distributed  with  varying 
statistical  character  over  range  and  bearing.  In  such  a  situation,  detectors  de¬ 
signed  under  a  Rayleigh  assumption  will  exhibit  increased  and  varying  false 
alarm  performance.  This  report  develops  a  technique  for  dealing  with  the 
unknown  and  varying  reverberation  statistics  by  using  a  non-Rayleigh  rever¬ 
beration  model  to  perform  a  statistical  normalization  of  the  background  rever¬ 
beration.  In  this  manner,  range-bearing  images  will  be  produced  that  have  a 
constant  background  reverberation  distribution  (i.e.,  Rayleigh).  The  technique 
is  evaluated  through  simulation,  with  particular  attention  to  the  false  alarm 
performance,  which  is  seen  to  depend  on  the  severity  of  the  non-Rayleighness 
of  the  reverberation  and  the  amount  of  data  used  to  estimate  the  parameters 
of  the  reverberation  model.  Application  to  real  data  has  shown  that  a  large 
degree  of  the  non- stationary  of  the  range-bearing  image  can  be  removed  by 
statistical  normalization. 

Keywords:  non-Rayleigh  o  reverberation  o  constant  false  alarm  rate  o  nor¬ 

malization 
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1 

Introduction 


The  detection  and  localization  of  submarines  using  active  sonar  is  often  hindered 
by  reverberation  from  the  ocean  surface,  volume  and  bottom.  In  shallow  water, 
reverberation  from  the  bottom  typically  dominates  the  other  types  of  reverbera¬ 
tion.  Normalization  algorithms  in  an  active  sonar  system  must  remove  the  non- 
stationarity  of  the  reverberation  prior  to  performing  detection.  These  algorithms 
traditionally  assume  the  reverberation  to  be  the  result  of  multiple  point  scatters, 
producing  a  Rayleigh  distributed  matched  filter  amplitude  with  a  time- varying  scale 
factor.  Thus,  the  only  effect  that  must  be  removed  is  the  time-varying  scale.  Many 
normalizers  have  been  designed  from  these  assumptions,  varying  primarily  in  their 
methods  for  dealing  with  corruption  of  the  auxiliary  data  used  to  estimate  the  mean 
power  level  [1,  2,  3,  4,  5].  However,  the  validity  of  the  multiple  point  scatterer  and 
resulting  Rayleigh  assumptions  will  not  hold  in  all  situations.  Many  researchers  (at 
a  variety  of  grazing  angles,  beamwidths,  frequencies,  transmit  signal  types,  and  bot¬ 
tom  conditions)  have  shown  that  bottom  reverberation  is  often  non-Rayleigh  (see  [6] 
for  a  comprehensive  bibliography).  Many  of  these  researchers  have  attempted  to  fit 
observed  reverberation  to  traditional  probability  distributions,  predominantly  the 
log-normal  distribution  [7,  8,  9,  10]  and  the  Weibull  distribution  [8,  9].  Others  have 
developed  their  own  models  including  one  based  on  a  Markov  driven  change  be¬ 
tween  two  different  types  of  sea-floor  patches  [11],  a  multiplicative  model  involving 
chi-squared  random  variables  [12]  and  a  mixture  or  multi-modal  Rayleigh  model  [13]. 
It  is  generally  accepted  that  non-Rayleigh  reverberation  is  caused  by  small  resolu¬ 
tion  cell  sizes  resulting  in  fewer  scatterers  so  that  the  central  limit  theorem  does  not 
hold  [14].  However,  it  has  also  been  observed  that  non- Rayleigh  reverberation  can 
arise  when  the  composition  of  the  sea-bed  varies  [11,  15]  or  when  Rayleigh  backscat- 
ter  is  affected  by  a  slow  modulation  as  in  radar  sea-surface  clutter  [16].  In  active 
sonar,  large  arrays  and  wide  bandwidth  transmit  signals  can  result  in  small  enough 
range-bearing  cells  to  produce  non-Rayleigh  reverberation.  In  low-frequency  sonar 
systems,  changes  in  the  bottom  composition  are  almost  inevitable  over  the  ranges 
of  interest. 

To  illustrate  the  non-Rayleigh  character  of  shallow  water  reverberation,  the  prob¬ 
ability  of  false  alarm  (Pfa)  measured  from  some  low  frequency  active  sonar  data 
obtained  during  SACLANT  Undersea  Research  Centre’s  SCARAB  trial  in  June, 
1997  is  compared  with  the  Rayleigh  distribution  in  Fig.  1.  If,  for  example,  a  de¬ 
tection  algorithm  was  designed  for  a  Pfa  of  10-3  assuming  Rayleigh  reverberation, 
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the  actual  Pfa  in  this  case  would  be  more  than  an  order  of  magnitude  greater.  This 
illustrates  the  effect  non-Rayleigh  reverberation  can  have  on  the  probability  of  false 
alarm  in  automatic  detection  algorithms.  When  the  non-Rayleigh  character  changes 
with  range  and  bearing,  whole  regions  in  a  normalized  range-beam  number  display 
can  appear  brighter  than  other  more  Rayleigh-like  regions.  The  range-bearing  dis¬ 
play  of  one  ping  of  SCARAB  data  shown  in  Fig.  2  illustrates  this  where  the  angular 
regions  to  either  side  of  broadside  to  the  array  (beam  27)  are  brighter  than  other 
areas.  This  non-stationarity  makes  target  detection  and  classification  difficult  both 
visually  and  by  automated  techniques  using  multiple  beam  information  such  as  im¬ 
age  processing. 


Figure  1  Probability  of  false  alarm  measured  from  SCARAB  data  compared  with 
the  Rayleigh  distribution. 

This  report  proposes  the  use  of  a  statistical  normalizer  to  convert  non-st  at  ionary 
non-Rayleigh  reverberation  back  to  the  Rayleigh  distribution.  This  is  accomplished 
by  approximating  the  distribution  of  the  reverberation  by  an  appropriate  model  and 
applying  a  nonlinearity  to  convert  the  reverberation  to  the  Rayleigh  distribution. 
This  concept,  along  with  the  Rayleigh  mixture  reverberation  model,  is  discussed  in 
Section  2.  The  Rayleigh  mixture  model  is  a  general  model  that  has  been  found  to 
adequately  represent  a  wide  variety  of  reverberation  [6,  15].  One  of  the  primary 
difficulties  in  implementing  such  a  statistical  normalizer  lies  in  obtaining  target- 
free  auxiliary  data.  These  data  are  required  to  estimate  parameters  describing  the 
distribution  of  the  diffuse  reverberation  and  must  not  be  corrupted  by  the  target 
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or  any  target-like  coherent  reverberation.  A  nonparametric  technique  for  removing 
strong  interferences  from  auxiliary  data  is  developed  in  Section  3.  A  new  method 
of  evaluating  the  false  alarm  performance  of  a  normalization  algorithm  is  developed 
in  Section  4  and  applied  to  both  the  statistical  normalizer  and  the  non-parametric 
interference  removal  algorithms  of  Sections  2  and  3.  The  performance  of  the  com¬ 
bined  statistical  normalizer /nonparametric  interference  pruning  algorithm  is  exam¬ 
ined  through  simulation  in  Section  5  and  then  applied  to  low  frequency  active  sonar 
data  obtained  from  SACLANT  Undersea  Research  Centre’s  SCARAB97  sea-trial  in 
Section  6. 


Matched  Filtered  Beam  Data 


Mean  Power  Level  Normalizer 


Range  (km) 

Figure  2  Range-beam  number  displays  of  the  matched  filtered  beam  data  (upper 
plot)  and  the  mean  level  normalizer  output  (lower  plot).  The  effect  of  non- Rayleigh 
reverberation  is  seen  in  the  normalizer  output  in  angular  regions  to  each  side  of 
broadside  to  the  array  (beam  27). 
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2 

Statistical  normalization 


2. 1  Concept 


Traditional  normalizers  assume  that  reverberation  is  Rayleigh  distributed  and  are 
designed  to  remove  the  changing  mean  power  level.  When  no  target  is  present, 
these  mean  power  level  normalizers  produce  a  homogeneous  or  identically  distributed 
background  with  a  known  probability  distribution.  Detection  is  then  straightforward 
to  perform  and  allows  exact  specification  of  the  Pja.  When  the  reverberation  is  non- 
Rayleigh,  a  nonlinearity  is  required  to  convert  it  to  being  Rayleigh  distributed. 


Suppose  the  observed  reverberation  data  Xi  have  probability  distribution  function 
(PDF)  f(x]Q)  and  cumulative  distribution  function  (CDF)  F(x]Q)  where  0  is  a 
parameter  vector.  Consider  the  random  variables  obtained  by  taking  the  value  of 
the  CDF  at  Xi, 


Ui  =  F(Xi-e).  (1) 

The  CDF  of  U{  is 

F\j  (u)  =  Pr  {U  <u} 

=  Pr  {F  (Xf,  0)  <  u} 

=  Pr  jAQ  <  F~l  (u;  0)} 

=  F  (f~x  (u;  ©) ;  ©) 

(2) 


for  u  E  (0, 1)  where  F~l  (u\  0)  is  the  functional  inverse  of  F  (x;  0),  which  is  assumed 
to  be  strictly  monotonically  increasing  over  the  region  (0,  oo)  to  provide  a  one-to- 
one  transformation  from  X%  G  (0,  oo)  to  Ui  G  (0, 1).  From  the  CDF  Ffj  (u)  —  u , 
it  is  clear  that  the  random  variables  Ui  are  uniformly  distributed  between  zero  and 
one;  that  is,  fu  (u)  —  1  for  u  e  (0, 1).  Similarly,  it  can  be  shown  that  the  random 
variables  taken  by  forming 

^  =  F^1  (Ui)  (3) 

will  have  CDF  Fr  ( y )  if  the  Ui  are  uniformly  distributed  between  zero  and  one, 
Fjj}  (u)  is  the  functional  inverse  of  the  proper  CDF  Fr  ( y ),  and  Fr  ( y )  is  strictly 
monotonically  increasing  over  (0,  oo). 
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In  this  manner  the  observed  reverberation  may  be  transformed  so  that  it  follows  the 
standard  Rayleigh  distribution  by  applying  the  non-linearity 

9  (*)  =  Fr1  (F  (*;  ©))  (4) 

for  x  >  0  where 

Fr  (y)  =  1  -  e~y2  for  y  >  0  (5) 

and 

Fr]  ( u )  =  yj-  log  (1  -  it)  for  u  e  (0, 1)  (6) 

are  respectively  the  Rayleigh  CDF  and  its  functional  inverse.  Now,  suppose  that 
the  reverberation  is  Rayleigh  with  power  A.  Then  the  CDF  is  simply 

F  (x;  A)  —  1  —  e~“  for  x  >  0  (7) 

and  the  nonlinearity  of  eq.  (4)  becomes 

9ix)  =  ^=  (8) 

which  is  identical  to  the  traditional  mean  level  normalizes  If  the  reverberation 
were  Weibull  distributed  with  CDF 


F  {x\  a,/3)  =  1  —  e  ax(3  for  x>0 


(9) 


then  the  nonlinearity  is 


r—  0 

g  (x)  =  yja  x  2 .  (10) 

To  illustrate  the  transformation  required  for  non-Rayleigh  reverberation,  the  non- 
linearities  required  to  convert  reverberation  that  is  Weibull,  K,  or  Rayleigh  mix¬ 
ture  distributed  are  shown  in  Fig.  3  where  the  de-emphasizing  of  the  larger  values, 
which  pulls  in  the  tails  of  the  reverberation,  is  clearly  seen.  The  reader  is  referred 
to  [6]  for  a  complete  description  of  these  non-Rayleigh  reverberation  models.  The 
parameters  for  the  non-Rayleigh  distributions  were  chosen  so  that  the  power  was 
one  and  the  kurtosis  equal  to  1.5  (the  kurtosis  of  a  Rayleigh  random  variable  is 
^4-^  ~  3  ~  0.2451).  The  Rayleigh  mixture  had  two  components  with  proportions 
0.9  and  0.1  and  powers  chosen  according  to  the  specified  kurtosis  and  total  power. 

From  a  statistical  perspective,  the  transformation  described  by  eq.  (4)  is  not  neces¬ 
sarily  the  most  advantageous  one.  Ideally,  the  log-likelihood  ratio  (LLR)  of  the  data 
would  be  used  to  implement  a  detector  that  maximizes  the  probability  of  detection 
while  constraining  the  probability  of  false  alarm.  However,  if  both  the  LLR  and 

SACLANTCEN  SR-303 


-  5  - 


Nonlinearity  Output 


SACLANTCEN  SR-303 


Figure  3  Example  tail  de- emphasizing  non-linearities  for  Weibull,  K,  and  Rayleigh 
mixture  distributions  with  a  kurtosis  value  of  1.5. 
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g  ( x )  are  monotonically  increasing,  then  either  may  be  used  if  the  detector  simply 
compares  each  individual  matched  filter  output  to  a  threshold  (i.e.,  no  combination 
of  the  data  over  time).  Conversely,  if  the  detector  combines  data  from  adjacent 
samples  in  time,  use  of  the  transformed  data  may  not  be  optimal. 

The  two  primary  difficulties  in  applying  the  nonlinearity  described  by  eq.  (4)  are 
first  choosing  the  appropriate  reverberation  model  F(x;0),  and  second,  estimating 
the  parameters  0.  The  Edgeworth  expansion,  log-normal,  Weibull,  K  and  Rayleigh 
mixture  distributions  and  their  parameter  estimation  are  discussed  in  [6]  along  with  a 
method  for  choosing  which  model  best  fits  the  observed  data.  One  of  the  conclusions 
of  [6]  was  that  the  Rayleigh  mixture  model  was  the  most  robust  in  fitting  different 
types  of  non- Rayleigh  reverberation,  a  result  also  supported  by  [15].  For  this  reason, 
the  Rayleigh  mixture  model  and  its  maximum  likelihood  parameter  estimation  are 
described  in  the  following  section. 


2.2  Rayleigh  mixture  model 

The  results  of  [7]— [13]  illustrate  the  variability  of  the  probability  distribution  of  the 
reverberation  amplitude.  A  single  model  that  is  general  enough  to  adequately  rep¬ 
resent  the  probability  distribution  of  any  encountered  reverberation  is  desired.  The 
lack  of  sufficient  and  accurate  information  about  sea- bottom  composition,  especially 
over  the  ranges  of  interest  to  low-frequency  active  sonar,  predicates  a  model  indepen¬ 
dent  of  environmental  parameters.  This  requirement  leads  to  considering  generalized 
noise  distributions  such  as  the  generalized  Gaussian,  generalized  Cauchy,  stable,  or 
mixture  distributions  [17].  Many  of  these  generalized  noise  models  describe  de¬ 
partures  from  a  nominal  distribution  which  then  becomes  a  special  case  (i.e.,  the 
generalized  Gaussian  model  has  as  special  cases  the  Laplace  and  standard  Gaus¬ 
sian  distributions);  however,  the  mixture  model  stands  out  in  its  flexibility.  The 
mixture  model  may  be  used  to  represent  an  impulsive  component  [17],  multi-modal 
data  [18],  and  can  accurately  represent  most  distributions  given  enough  components 
in  the  mixture. 

Under  the  assumption  that  the  majority  of  active  sonar  reverberation  will  result  in 
a  near  Rayleigh  distribution,  a  mixture  of  Rayleigh  random  variables  with  different 
variances  may  be  general  enough  to  adequately  describe  the  reverberation  PDF  with 
a  minimal  parameterization.  This  is  somewhat  supported  by  the  results  of  Stewart 
et  al.  [13]  where  good  agreement  was  found  between  a  Rayleigh  mixture  and  the  ob¬ 
served  reverberation  data  and  by  those  of  Trunk  [19]  where  two  component  Gaussian 
mixtures  were  used  to  model  the  quadrature  and  in-phase  parts  of  the  matched  filter 
output  for  radar  clutter.  Stronger  support  may  be  found  in  Crowther’s  [11]  model 
of  backscatter  from  a  patchy  seabed  which  simplifies  to  a  two  component  Rayleigh 
mixture  if  the  resolution  cell  size  is  small  compared  with  the  patch  size.  A  further 
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and  more  mundane  justification  for  considering  a  Rayleigh  mixture  model  lies  in 
the  ability  to  obtain  the  maximum  likelihood  estimates  of  the  mixture  parameters 
through  the  expectation-maximization  (EM)  algorithm  [20,  21]. 

The  Rayleigh  PDF  has  form 

/  (®;  A)  =  ye‘T  (11) 

for  x  >  0  where  A  is  the  power  (i.e.,  A  =  E  [x2]).  A  mixture  of  p  Rayleigh  random 
variables  has  the  PDF 

f(x;G)  =  /(^Ai)  (12) 

i= 1 

for  x  >  0  where  0  =  •  •  •  %  Ai  •  •  •  Ap],  7 r*  is  the  mixture  proportion  and  A * 

is  the  power  of  the  component.  The  corresponding  CDF  is 

p 

F  (x;  0)  = 

i=l 

P  x2 

=  ^7^6  Xi  (13) 

i=l 

for  x  >  0. 

Substitution  of  eq.  (13)  into  the  nonlinearity  of  eq.  (4)  results  in 

g{x)  =  -logj^TT.e  |  .  (14) 

It  should  also  be  noted  that  the  standard  Rayleigh  distribution  results  if  there  is 
only  one  component  in  the  Rayleigh  mixture  or  if  all  components  have  the  same 
power. 

Care  should  be  taken  in  the  numerical  evaluation  of  eq.  (14),  particularly  for  large 
values  of  x  where  finite  precision  arithmetic  may  return  a  value  of  zero  for  the 
summation  even  though  the  summand  is  always  positive  (for  7r,  ^  0).  A  solution  is 
to  factor  out  the  largest  term  in  the  summation,  say  7r0exp  which  simplifies 

further  once  the  logarithm  is  applied.  The  resulting  summation  is  then  one  plus 
p  —  1  terms  between  zero  and  one  and  may  be  evaluated  without  fear  of  significant 
numerical  error. 

2.2.1  Target  models 

Another  distinct  advantage  of  the  Rayleigh  mixture  model  is  it’s  ability  to  incor¬ 
porate  the  standard  fluctuating  and  non-fluctuating  target  models.  The  Rayleigh 
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mixture  model  may  be  developed  from  a  mixture  of  zero-mean  complex  Gaussian 
random  variables  for  the  quadrature  and  in-phase  components  of  the  matched  filter 
output.  This  differs  slightly  from  the  work  of  Trunk  [19]  in  that  this  model  forces 
both  the  quadrature  and  in-phase  components  to  come  from  the  same  mixture  com¬ 
ponent  and  results  in  a  Rayleigh  mixture  for  the  magnitude. 


The  non-fluctuating  target  model  is  represented  by  a  non-zero  mean  component  in 
the  complex  matched  filter  output— essentially  adding  a  constant  to  the  aforemen¬ 
tioned  mixture  of  zero-mean  complex  Gaussian  random  variables.  The  resulting 
magnitude  distribution  is  a  mixture  of  Rician  random  variables 


p 

/nf  0;Q>m)  =  h 


i—1 


(15) 


where  js  a  signal-to- reverberation  power  ratio  (SRR), 


Onr 

fi  (x;\ ,6)  =  —  e 


(16) 


is  the  Rician  PDF  with  scale  A  and  SRR  parameter  <5,  and  Iq  ( x )  is  the  zero-order 
modified  Bessel  function. 


The  fluctuating  target  model  is  represented  by  an  additive  zero-mean  Gaussian 
component  to  the  complex  matched  filter  output.  This  results  in  a  mixture  of 
Rayleigh  random  variables  with  power  A*  =  A^  +  Ao  where  Ao  is  the  variance  of  the 
target  component, 

/a  (®;  ©>  Ao)  =  y /  (^;  Ai  +  Ao) .  (17) 

i= 1 

The  proof  of  these  results  is  sketched  in  Annex  A. 


2.2.2  Estimation  of  the  Rayleigh-mixture  parameters 

The  Rayleigh  PDF  falls  under  the  exponential  family  of  distributions  [22].  Thus, 
the  EM-algorithm  may  be  applied  directly  from  the  results  of  [20]  or  [21]  to  obtain 
the  maximum  likelihood  estimates  of  the  mixture  parameters.  For  the  interested 
reader,  this  development  is  provided  here  in  Annex  B. 

The  EM-algorithm  for  obtaining  maximum  likelihood  estimates  of  the  mixture  pa¬ 
rameters  from  the  observed  data  (xi,  . ..,  x„)  is  summarized  as  follows: 

1.  Initialize  ©  -  For  example,  set  7r,  =  |  and  form  the  A*  as  the  sample  power 
of  subsets  of  the  observed  data. 
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2.  Expectation  step  -  Form  the  intermediate  variables 


4 

A  ?Le-t 

J  A, 

(18) 

and 

‘<>5 

W 

li 

*«*■» 

(19) 

for  i  =  1,  . . . ,  p  and  j  =  1,  . . . ,  n. 

3.  Maximization  step  -  Form  new  estimates  of  0 

3= 1 

(20) 

and 

\  ' J— 1  3  LiJ 

1  ~  yn  ,  w •  •  ‘ 

rr  i,j 

(21) 

4.  Convergence  -  If  the  normalized  change  between  consecutive  estimates  of  the 
unknown  parameters  is  greater  than  some  tolerable  value,  then  go  to  step  2. 


As  the  EM  algorithm  converges  to  local  and  not  global  maxima,  it  may  be  ap¬ 
propriate  to  run  the  iteration  from  several  initial  starting  points.  A  reasonable 
initialization  for  the  mixture  proportions  is  to  set  them  to  be  equally  likely,  ki  = 
Reasonable  initial  choices  for  the  power  components  are  either  (i)  the  sample  pow¬ 
ers  formed  from  p  disjoint  subsets  of  the  n  data  samples,  or  (ii)  the  sample  powers 
formed  from  p  subsets  of  the  n  ordered  data  samples  so  that  A i  <  A^+i. 

Choice  of  the  model  order  poses  a  slightly  more  difficult  problem.  It  is  possible  to 
apply  Akaike’s  information  criterion  [23]  or  Rissanen’s  minimum  descriptive  length 
criterion  [24]  to  choose  the  model  order.  Other  equally  justifiable  methods  include 
evaluation  of  the  Kolmogorov- Smirnov  statistic  or  consideration  of  the  mixture  com¬ 
ponent  proportions  and  powers  (i.e.,  if  an  additional  component  results  in  a  near 
zero  proportion  or  a  power  identical  to  a  component  already  present,  the  lower  order 
model  is  likely  adequate).  Experience  has  shown  that  three  or  fewer  components 
are  usually  sufficient  to  represent  low-frequency  active  sonar  reverberation. 

In  applying  this  method  to  real  data,  the  time  required  for  convergence  becomes 
an  issue  owing  to  the  massive  number  of  times  the  parameter  estimation  must  be 
performed  (essentially  once  for  every  range-bearing  cell).  As  it  is  expected  that  the 
non-Rayleigh  character  of  the  data  will  not  change  too  rapidly  in  range  or  bearing, 
the  parameter  estimation  may  not  need  to  be  performed  for  every  range  sample  and 
the  time  to  convergence  may  be  reduce  by  initializing  the  algorithm  with  the  most 
recent  set  of  parameter  estimates. 
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3 

Removal  of  interference  from  auxiliary  data 


In  constant  false  alarm  rate  (CFAR)  detection,  masking  of  a  target  occurs  when  the 
estimate  of  the  reverberation  power  formed  from  the  auxiliary  data  is  unnaturally 
high,  resulting  in  a  reduced  probability  of  detection.  This  occurs  particularly  when  a 
target  or  target-like  interference  is  present  in  the  auxiliary  data  and  results  in  shadow 
regions  around  the  interference.  Normalizers  based  on  order  statistics  (OS)  [2]  are 
traditionally  used  to  alleviate  problems  such  as  masking  or  its  counterpart,  false 
alarm  rate  inflation.  These  normalizers  use  a  subset  of  the  ordered  auxiliary  data 
to  form  an  estimate  of  the  reverberation  power,  essentially  discarding  some  of  the 
largest  and  smallest  values.  However,  they  do  not  utilize  any  information  about  the 
location  of  the  largest  or  smallest  values  within  the  auxiliary  data  window.  Lops 
and  Willett  [4]  have  proposed  a  CFAR  detector  that  does  exploit  position  informa¬ 
tion  using  the  LZ-filters  of  Palmieri  and  Boncelet  [25].  Their  LZ-CFAR  detector  is  a 
mean  power  level  normalizer  with  robust  performance  when  the  auxiliary  data  are 
corrupted  by  interference  (target-like  data)  and  clutter  (regions  with  higher  inten¬ 
sity).  However,  as  previously  noted,  when  dealing  with  non-Rayleigh  reverberation, 
normalization  requires  more  than  simple  estimation  of  the  mean  in  the  presence  of 
interferences  in  the  auxiliary  data.  An  uncorrupted  sample  of  the  reverberation  is 
required  to  determine  the  departure  of  the  distribution  from  Rayleigh. 

If  strong  range-spread  target-like  interferences  are  present  in  the  auxiliary  data, 
removal  of  some  of  the  largest  values  would  successfully  produce  a  homogeneous 
(i.e.,  uncorrupted)  sample  of  reverberation-only  data.  However,  when  the  auxiliary 
data  are  homogeneous,  such  a  process  would  remove  the  largest  reverberation  values. 
Removal  of  these  extreme  values  eliminates  evidence  of  heavy  tails  from  the  auxiliary 
data  and  would  result  in  a  flawed  estimation  of  the  PDF  of  the  auxiliary  data.  Thus, 
an  alternative  method  is  required  to  remove  range-spread  target-like  interferences 
when  they  are  present  and  leave  the  auxiliary  data  unmodified  otherwise.  It  is 
possible  to  prune  these  interferences  from  the  auxiliary  data  by  exploiting  their 
spatial  compactness  (in  range).  This  is  also  the  premise  of  the  LZ-CFAR  detector 
of  Lops  and  Willett  [4]  where  an  LZ-filter  is  used  to  estimate  the  mean  power  level. 
Here,  however,  a  decision  is  made  first  about  the  presence  of  an  interference  in  the 
auxiliary  data  followed  by  removal  of  the  affected  data  only  if  it  is  decided  that  an 
interference  is  present.  The  resulting  pruned  auxiliary  data  is  then  used  to  estimate 
the  pertinent  reverberation  model  parameters. 
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As  the  PDF  of  the  reverberation  will  not  be  known  exactly  and  is  expected  to  change 
with  range  and  bearing,  a  nonparametric  technique  is  required  to  choose  between 
the  interference-present  hypothesis  (Hx)  and  the  null  hypothesis  (Ho)  of  indepen¬ 
dent  and  identically  distributed  reverberation  data.  Let  the  set  {KX,K2,  ■ •  ■ ,  Kn} 
represent  the  integer-valued  indices  to  the  ordered  complete  auxiliary  data.  If 
{Xi,X2,...,Xn}  are  the  n  samples  of  auxiliary  data,  then  the  ordered  samples 

are 

X(1)  <  A(2)  <•••  <  X{n)  (22) 

and  if  X{  =  X(n_J+1)  then  Kj  =  i.  For  simplicity,  the  indices  have  been  numbered 
contrary  to  standard  notation  for  ordered  statistics  so  that  Kx  represents  the  position 
or  index  of  the  largest  data  sample,  K2  the  second  largest,  down  to  Kn  representing 
the  position  of  the  smallest  sample.  If  the  auxiliary  data  are  statistically  independent 
and  identically  distributed,  then  the  probability  of  the  largest  value  occurring  in  the 
first  position,  or  second,  or  nth  are  all  equally  that  is,  the  probability  mass 
function  (PMF)  of  K\  is 

V  (&i)  =  -  for  ki  =  1, . . . ,  n.  (23) 

n 

Similar  arguments  lead  to  the  marginal  PMF  of  any  one  index 

p  (kj)  =  t  for  kj  =  1,  •  •  • ,  n  (24) 


and  the  joint  PMF  of  any  pair  of  indices 


p(h,kj)  - 


for  kt  =  1,  •  •  • ,  n  and  kj  =  1,  . . . ,  n  excepting  cases  where  h  =  kj.  Extension  to  the 
joint  PMF  of  all  the  indices  results  in 

p(kh...,kn)  =  ^  (26) 


where  each  index  varies  between  1  and  n  and  no  two  indices  can  have  the  same 
value.  This  result  implies  that  any  test  statistic  formed  from  the  indices  to  the 
ordered  statistics  will  be  nonparametric  under  the  null  hypothesis. 


Formulation  of  an  optimal  test  is  not  possible  due  to  the  vague  description  of  the 
interference-present  hypothesis  (Hx)  and  the  difficulty  in  analyzing  the  resulting 
probability  mass  function  (PMF)  of  the  indices  to  the  ordered  statistics.  However, 
if  an  interference  is  spatially  compact  and  stronger  than  the  underlying  noise  dis¬ 
tribution,  it  is  likely  that  the  indices  of  the  largest  order  statistics  will  indicate  the 
position  of  such  an  interference  within  the  auxiliary  data.  This  information,  accom¬ 
panied  by  the  fact  that  all  the  permutations  of  the  indices  to  the  order  statistics  are 
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equally  likely  under  the  null  hypothesis,  leads  to  considering  the  distance  between 
the  indices  pointing  to  the  largest  values  of  the  auxiliary  data.  Suppose  only  the 
indices  to  the  two  largest  samples  are  considered.  The  test  statistic 

T2  =  (K2  -  Kif  (27) 

may  provide  adequate  distinction  between  Ho  and  Hi  if  both  of  the  two  largest 
samples  are  due  to  the  interference.  With  heavy  tailed  reverberation,  target-like 
bottom  features  of  unknown  size,  and  spreading  of  the  target  return  due  to  the 
multipath  of  shallow  water  and  reflection  off  of  a  large  target,  this  may  not  always 
be  the  case.  Thus,  it  is  prudent  to  also  consider  the  index  to  the  third  largest  sample 
and  form  a  similar  test  statistic, 

T3  =  ( if  2  —  -Ki )2  +  (K%  —  K%)  +  —  K\) 

=  T2  +  {Kz-K2)2  +  {Kz-Kif ,  (28) 

which  will  perhaps  perform  well  if  at  least  two  of  the  three  largest  samples  are  due 
to  the  interference.  Extending  this  type  of  test  statistic  to  one  considering  the  m 
largest  samples  results  in 

m  i— 1 

Tm  =  EE(^-^)2 
2  1 

m— 1 

=  Tm- 1  +  ^  (Km  —  Ki)  .  (29) 

i=l 

If  a  strong  and  spatially  compact  interference  exists  in  the  auxiliary  data,  the  pruning 
statistics  Tm  should  be  smaller  than  when  the  auxiliary  data  are  uncorrupted.  Thus, 
interference  presence  may  be  tested  by  choosing 

H0  :  Interference  absent  if  Tm>  hm,  or 

Hi  :  Interference  present  if  Tm  <  hrn 

where  hm  is  a  threshold  that  is  ideally  chosen  according  to  a  desired  probability  of 
falsely  declaring  an  interference  present  in  the  auxiliary  data  (i.e.,  Type  I  error). 
As  the  spatial  extent  of  the  interference  is  not  known  a  priori,  it  may  be  difficult 
to  choose  an  appropriate  value  for  m.  As  an  alternative,  each  value  of  Tm  may  be 
considered  up  to  some  maximum  number  mmax  which  should  be  much  smaller  than 
n  (Tn  is  a  constant)  and  greater  than  the  maximum  spatial  extent  of  an  interference. 
Unfortunately,  this  greatly  complicates  analysis  of  the  Type  I  error  probability  and 
thus  the  thresholds  hm.  In  fact,  computation  of  the  PMF  of  Tm  under  Ho  is  reason¬ 
ably  feasible  analytically  only  for  m  =  2  and  numerically  so  only  for  small  m  and  n. 
Thus,  it  is  not  possible  to  exactly  determine  thresholds  hm  based  on  the  probability 
that  Tm  <  hm  under  H0  for  anything  but  small  m. 

As  the  consequence  of  falsely  declaring  an  interference  present  in  the  auxiliary  data 
is  the  removal  of  the  noise  outliers,  exact  knowledge  of  the  Type  I  error  probability 

SACLANTCEN  SR-303 


-13- 


SACLANTCEN  SR-303 


is  not  crucial.  Thus,  approximate  thresholds  may  be  obtained  by  assuming  a  Gaus¬ 
sian  form  for  the  PMF  of  Tm.  The  Gaussian  PDF  is  completely  determined  by  its 
mean  and  variance,  which  in  this  case  entails  a  fairly  involved  derivation  as  seen  in 
Annex  C.  The  results  are  as  follows, 


E  [Tm] 


'm  (m  —  1)' 

~n  (n  +  1)' 

2 

6 

J  \n{n+  1) 

1 L  6 


(30) 


and 


Var  [Tm] 


where 


l  {2c\  —  8c2  +  6C3)  +  l0  (4c3  +  4c4  —  8C5) 


+  [l  (Z  —  1)  —  Jo]  (ci  —  4c2  +  3C3)  —  Z2 


n2  (n  +  1)^ 


36 


(31) 


1  = 


m  (m  —  1) 


(32) 


Zo  = 


6(T)  rn>  4 
0  m  <  4 


(33) 


and  C1-C5  are  as  described  in  Annex  C.  The  threshold  hm  is  then  chosen  according 
to  a  desired  probability  of  Type  I  error,  a,  under  a  Gaussian  assumption  for  Tm, 

hm  —  E  [. Tm \  —  za  ^/Var  [T^],  (34) 

where  is  the  a-quantile  of  a  Gaussian  random  variable  (Z)  with  zero  mean  and 
unit  variance, 


a  —  Pr  {Z  >  Zq)  =  Pr  {Z  <  - z Q}  . 


(35) 


3.1  Choosing  which  cells  to  remove 

Once  the  decision  has  been  made  that  a  spatially  compact  interference  is  corrupting 
the  auxiliary  data,  the  affected  cells  must  be  removed.  Without  exact  knowledge 
about  the  range  extent  of  the  interference  it  is  difficult  to  know  how  many  and  which 
cells  to  remove.  The  following  nonparametric  technique  exploits  the  expected  spa¬ 
tial  compactness  of  the  interference,  additionally  allowing  the  removal  of  a  variable 
number  of  cells. 


1.  Isolate  the  indices  pointing  to  the  m  largest  values  in  the  auxiliary  data  sample. 
Choose  m  according  to  the  maximum  number  of  possible  interference  samples 
(perhaps  rnmax). 
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2.  Form  a  histogram  using  equally  spaced  bins  between  the  smallest  and  largest 
indices  of  the  interference  subset.  The  number  of  bins  should  be  less  than  or 
equal  to  m. 

3.  Discard  any  cells  within  a  reasonable  distance  (say  half  to  two-thirds  of  the 
maximum  number  of  interference  cells)  of  the  cell  yielding  the  largest  frequency 
from  the  histogram. 

3.2  Discussion 

The  non-parametric  interference  pruning  (NPIP)  algorithm  described  in  the  previ¬ 
ous  sections  was  developed  to  deal  with  the  existence  of  a  strong  spatially  compact 
interference  with  unknown  range  extent  that  corrupts  possibly  non-Rayleigh  aux¬ 
iliary  data.  The  objective  is  to  produce  an  interference-free  sample  of  the  rever¬ 
beration  data  from  which  reverberation  model  parameters  may  be  estimated.  The 
requirement  for  the  procedure  to  deal  with  non-Rayleigh  reverberation  with  un¬ 
known  distribution  predicates  a  non-parametric  technique.  The  strong  interference 
assumption  is  arrived  at  by  postulating  that  weak  interferences  will  not  greatly  alter 
the  reverberation  parameter  estimation  and  additionally  helps  to  justify  the  non- 
parametric  form  of  the  algorithm.  The  unknown  range  extent  of  the  interference  is 
required  to  account  for  interferences  (i.e.,  bottom  features  producing  coherent  rather 
than  diffuse  reverberation)  of  unknown  dimension  and  also  for  the  spreading  that  an 
active  sonar  echo  incurs  from  propagation  to  and  from  the  target  and  reflection  off 
of  the  target.  This  particular  assumption  is  certainly  valid  for  current  tactical  sonar 
bandwidths  and  array  sizes;  however,  it  may  need  to  be  reconsidered  if  no  spreading 
is  expected  (i.e.,  if  the  target  and  interference  are  smaller  than  the  resolution  cell 
size). 

The  NPIP  algorithm  suffers  from  one  glaring  weakness — the  lack  of  consideration  of 
multiple  targets  or  interferences  within  the  auxiliary  data  set.  Accounting  for  more 
than  one  target  or  interference  greatly  complicates  the  problem.  Possible  solutions 
include  applying  the  NPIP  algorithm  to  windows  small  enough  that  the  existence 
of  multiple  interferences  is  not  likely  (i.e.,  apply  it  individually  to  the  leading  and 
lagging  windows  of  auxiliary  data)  or  performing  a  pre-detection  screening  of  the 
data  using  a  non-parametric  test  like  the  Page  test  described  in  [26]  or  by  passing 
a  small-window-NPIP  over  the  time  series  and  tagging  all  strong  spatially  compact 
interferences.  These  strong  deviations  from  the  background  would  then  be  excluded 
from  use  in  the  normalization  process. 
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4 

Instantaneous  probability  of  false  alarm 


The  primary  goal  of  a  normalization  scheme  is  to  provide  a  constant  false  alarm 
probability  in  the  presence  of  a  specific  type  of  change  in  the  null  distribution, 
which,  in  this  case,  is  the  reverberation  distribution.  For  example,  the  mean  power 
level  normalizers  provide  constant  Pfa  for  any  reverberation  power  level.  There 
is,  however,  a  price  to  be  paid  for  this  invariance.  The  variability  of  parameter 
estimation  essentially  introduces  a  variability  on  the  threshold  used  for  detection, 
changing  both  the  detection  and  false  alarm  performance. 

In  this  section,  a  new  method  for  evaluating  and  comparing  the  false  alarm  perfor¬ 
mance  of  normalizers  is  proposed.  Traditionally,  normalizers  are  evaluated  by  their 
loss  in  detection  performance  when  the  detector  thresholds  are  chosen  to  provide 
the  same  Pfa.  However,  when  the  auxiliary  data  are  corrupted  by  interferences  or 
if  the  thresholds  are  chosen  based  on  incorrect  assumptions  about  the  reverberation 
statistics,  the  false  alarm  performance  changes  and  should  be  considered  separately 
and  in  addition  to  P&.  The  method  proposed  in  this  section  specifically  considers 
the  variability  on  the  false  alarm  performance  introduced  by  the  estimation  of  re¬ 
verberation  parameters  and  provides  an  idea  of  what  the  detection  performance  will 
be.  Detection  performance  is  considered  more  fully  in  Section  5. 


4.1  Mean  power  level  normalizer 


The  effect  of  the  variability  of  parameter  estimation  on  the  probability  of  false  alarm 
for  a  mean  level  normalizer  may  be  described  by 


Pfa 


Et 


(36) 


where  h  is  a  threshold  chosen  to  provide  a  specified  Pfa,  the  probability  is  taken 
under  the  null  or  reverberation-only  hypothesis  and  the  expectation  is  taken  over 
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the  reverberation  power  level  estimate  A.  The  quantity 

Pfa  (A,  h)  =  Pro  jx  >  hyfl  |  A  J 

=  1  —  F  (h\fl ;  a)  (37) 

where  F  (x]  A)  is  the  CDF  of  X  under  the  null  hypothesis  may  be  thought  of  as 
the  conditional  Pfa  in  a  manner  analogous  to  the  conditional  signal~to-noise  ratio 
of  Reed,  Mallett  and  Brennan  [27]  or  as  the  instantaneous  Pfa  in  the  sense  that  it 
is  the  Pfa  given  the  current  estimate  of  the  reverberation  power.  The  instantaneous 
Pfa  provides  a  measure  of  the  quality  of  the  parameter  estimation  and  an  indication 
of  the  effect  on  detection  performance,  particularly  in  the  presence  of  interferences 
in  the  auxiliary  data  window.  As  an  example,  consider  Fig.  4  where  the  PDF  of 
logio  (-P/a  (A, /i))  is  shown  for  a  cell- averaging  constant  false  alarm  rate  detector 
(CA-CFAR),  an  order  statistic  constant  false  alarm  rate  detector  (OS-CFAR),  and 
a  non-parametric  interference  pruning  constant  false  alarm  rate  detector  (NPIP- 
CFAR).  The  OS-CFAR  detector  used  the  82nc*  ordered  statistic  of  the  100  samples 
of  Rayleigh  distributed  auxiliary  data.  The  NPIP-CFAR  detector  removed  auxiliary 
data  potentially  corrupted  by  a  target  or  target-like  interference  by  assuming  a 
maximum  of  mmax  =  10  target  samples  and  using  thresholds  generated  according 
to  eq.  (34)  with  a  =  0.1.  The  remaining  data  were  then  used  in  the  manner  of  a 
CA-CFAR  detector  to  estimate  the  reverberation  power  level.  All  of  the  detectors 
have  their  thresholds  set  so  that  Pja  =  10“4,  for  which  the  NPIP-CFAR  required 
simulation  analysis.  As  expected,  the  PDFs  of  the  instantaneous  Pfa  for  the  OS- 
CFAR  and  NPIP-CFAR  have  broader  tails  than  those  of  the  CA-CFAR  detector 
indicating  slightly  worse  performance.  Note  that  the  actual  Pfa  is  the  expected 
value  of  10  raised  to  the  random  variable  represented  by  the  PDFs  shown  in  the 
figures. 

If  parameter  estimation  were  perfect,  the  PDF  of  the  instantaneous  Pfa  would  have 
an  impulse  function  at  the  actual  Pfa  and  the  corresponding  CDF  would  look  like 
a  rising  step  function  with  the  transition  at  the  Pja.  From  this  ideal,  one  could 
use  the  slope  of  the  CDF  of  the  instantaneous  Pfa  at  the  desired  Pfa  as  a  measure 
of  effectiveness.  The  slope  of  the  CDF  is,  of  course,  the  PDF,  and  indicates  the 
steepness  of  the  CDF  in  the  region  of  the  Pfa — a  steeper  CDF  at  the  mean  implies 
less  variance  in  the  PDF  and  thus  more  precise  false  alarm  performance.  This  was, 
in  fact,  how  the  82n^  order  statistic  was  chosen  in  the  previous  example,  and  it  is 
clearly  seen  in  Fig.  4  that  the  CA-CFAR  has  the  largest  PDF  at  the  desired  P/a. 

Introduction  of  a  non-fluctuating  (Swerling  Type  0  resulting  in  Rician  distributed 
data)  target-like  interference  into  5  of  the  100  samples  of  auxiliary  data  shifts  all 
of  the  instantaneous  Pfa  PDFs  to  the  left,  as  seen  in  Fig.  5,  indicating  a  reduction 
in  detection  performance  (target  masking).  The  interference  had  an  interference- 
to-reverberation  power  ratio  (IRR)  of  10  dB  and  was  randomly  positioned  within 
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the  auxiliary  data  for  each  of  the  10,000  simulation  trials.  Clearly  the  PDF  of  the 
instantaneous  Pfa  for  the  NPIP-CFAR  detector  is  affected  the  least  of  the  three  de¬ 
tectors,  remaining  largely  unchanged  from  the  case  where  no  interference  is  present. 
As  a  smaller  instantaneous  Pfa  implies  a  larger  effective  threshold  and  thus  a  smaller 
probability  of  detection,  clearly  the  NPIP-CFAR  detector  provides  the  best  perfor¬ 
mance  in  the  presence  of  such  an  interference  in  the  auxiliary  data. 


4.2  Statistical  normal izer 

The  instantaneous  Pja  may  be  similarly  developed  from  the  definition  of  the  Pfa  for 
the  statistical  normalizer  of  eq.  (4) 

Pfa  =  Pro  {</(*;©)  >h} 

=  Pro  f  X  >  g~l  (h\  ©)  } 

=  Eg,  [Pr0  {-X-  >  g~l  [h\  ©)  |©  }j 
-  Eg  {Pfa  (©,/*)]  (38) 

where  the  dependence  of  g  on  0  has  been  indicated  explicitly  and  g  1  is  the  func¬ 
tional  inverse  of  g.  The  instantaneous  Pfa  is  then 

Pfa  (©,  hj  =  Pro  { X  >  g-1  (h;  ©)  |©  } 

=  i  -  f  (g~x  (fy  ©) ;  ©)  (39) 

where  F  (x;  0)  is  the  CDF  of  X  under  the  null  hypothesis. 

The  instantaneous  Pfa  provides  a  complete  picture  of  the  false  alarm  performance, 
with  the  average  value  representing  the  standard  probability  of  a  false  alarm.  In 
the  presence  of  non-Rayleigh  reverberation,  the  consequences  of  assuming  Rayleigh 
reverberation  or  attempting  to  account  for  the  departure  from  Rayleigh  become 
evident  in  the  PDF  of  the  instantaneous  Pfa  as  do  the  effects  of  target  masking  and 
false  alarm  rate  inflation. 


4.3  Statistical  normalizer:  Weibull  reverberation 


If  X  is  Weibull  distributed  with  parameters  0  =  (a,  (3)  and  g  is  as  in  eq.  (10)  with 
the  parameter  estimates  0  =  ,  the  instantaneous  Pfa  is 


Pfa 


(40) 
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Figure  4  Probability  density  function  of  log10  of  instantaneous  Pfa  for  CA-CFAR , 
anrf  NPIP-CFAR  with  100  samples  of  Rayleigh  distributed  auxiliary  data . 
The  thresholds  for  each  detector  were  chosen  so  that  the  average  Pja  —  10“4. 


Figure  5  Same  as  the  above  figure  except  with  a  10  dB  IRR  non-fluctuating  target 
corrupting  5  samples  of  the  auxiliary  data .  The  actual  Pja  is  indicated  on  each  curve 
by  an  *. 
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If  the  parameters  are  estimated  perfectly  (i.e.,  0  =  0),  the  instantaneous  P/a  be¬ 
comes  constant  at  Pfa  ( h )  =  exp  {— /i2},  exactly  that  for  the  Rayleigh  distribution. 

In  [6]  the  maximum  likelihood  (ML)  estimators  for  the  Weibull  parameters  had  less 
error  between  the  parameter  estimate  generated  CDF  and  the  actual  CDF  but  did 
not  fit  the  sample  CDF  as  well  as  the  method  of  moments  (MM)  estimators  for 
moderate  kurtosis  values.  The  instantaneous  Pfa  provides  a  means  for  compar¬ 
ing  the  performance  of  these  estimators  in  a  statistical  normalizer  with  the  CA- 
CFAR  normalizer  using  a  measure  relevant  to  detector  performance.  The  PDFs  of 
logio  ( Pfa  (©5^))  of  eq.  (40)  for  the  statistical  normalizer  using  the  ML  and  MM 
parameter  estimators  and  for  the  CA-CFAR  detector  are  shown  in  Fig.  6.  The  PDFs 
are  estimated  from  10,000  simulation  trials  with  Weibull  distributed  reverberation 
when  the  power  and  kurtosis  are  both  1  (a  =  1.1003  and  (3  =  1.6147).  The  thresh¬ 
olds  for  all  of  the  detectors  were  chosen  assuming  Rayleigh  distributed  reverberation 
(i.e.,  h  =  y/—  logP/a).  The  PDFs  of  the  instantaneous  Pja  for  the  statistical  nor¬ 
malizes  using  the  ML  and  MM  parameter  estimates  are  nearly  indistinguishable 
and  clearly  provide  better  performance  than  the  CA-CFAR  detector  which  blindly 
assumes  the  data  to  be  Rayleigh.  Improved  false  alarm  performance  is  illustrated 
by  the  decreasing  variance  and  trend  of  the  mean  toward  the  desired  Pfa  as  more 
data  are  used  to  estimate  the  Weibull  parameters. 

The  actual  Pfa  and  the  standard  deviation  of  the  instantaneous  Pfa  are  shown  in 
Table  1  where  it  is  seen  that  the  MM  parameter  estimates  provide  slightly,  though 
not  significantly,  better  performance  for  these  particular  a  and  (3.  Increasing  the 
kurtosis  resulted  in  identical  performance  for  the  detector  using  the  ML  parameter 
estimates  and  slightly  worse  performance  for  the  detector  using  the  MM  parameter 
estimates.  The  equivalence  of  the  performance  over  varying  a  and  (3  for  the  detector 
using  the  ML  parameter  estimates  indicates  that  there  may  be  an  invariance  to  these 
parameters,  providing  a  constant  false  alarm  rate  over  the  complete  Weibull  family 
of  distributions.  This  invariance  does  in  fact  exist  when  the  ML  parameter  estimates 
are  used  in  either  the  nonlinearity  of  eq.  (10)  or  the  instantaneous  Pja  of  eq.  (40)  as 
shown  in  Annex  D.  The  importance  of  this  result  is  that  if  the  reverberation  is  truly 
Weibull  distributed,  then  the  false  alarm  performance  of  the  statistical  normalizer 
using  the  ML  parameter  estimates  and  the  non-linearity  of  eq.  (10)  is  constant 
regardless  of  the  values  of  a  and  /?.  Thus,  the  threshold  h  may  be  chosen  solely  as 
a  function  of  the  desired  Pfa  and  the  amount  of  auxiliary  data  used  to  estimate  the 
parameters  of  the  Weibull  distribution.  As  previously  indicated,  the  detector  using 
the  MM  parameter  estimates  does  not  enjoy  this  freedom.  Unfortunately,  the  lack  of 
a  closed  form  solution  for  the  ML  estimates  of  a  and  (3  precludes  an  exact  analytical 
analysis  of  the  false  alarm  performance,  thus  forcing  the  use  of  simulation.  The 
invariance  of  the  nonlinearity  of  eq.  (10)  to  changes  in  a  and  (3  provides  a  strong 
argument  for  using  the  ML  parameter  estimators  over  the  MM  estimators,  even 
though  they  may  be  slightly  more  computationally  intensive.  Goldstein  [28]  has 
also  proposed  a  detector  having  invariant  false  alarm  performance  over  the  family  of 
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Figure  6  Probability  density  function  of  log10  of  instantaneous  Pfa  for  statistical 
normalizers  using  ML  and  MM  parameter  estimates,  which  are  nearly  identical, 
and  for  CA-CFAR  detector  with  100,  500  and  1000  samples  of  Weibull  distributed 
auxiliary  data .  The  actual  Pja  for  each  detector  is  indicated  by  the  *. 
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Weibull  distributions,  though  it  does  not  have  the  maximum  likelihood  parameter 
estimation  argument  that  supports  the  statistical  normalizer. 


Table  1  Mean  (Pfa)  and  standard  deviation  of  instantaneous  Pfa  for  ML  and  MM 
parameter  estimates  of  Weibull  distributed  reverberation  using  auxiliary  data  sets  of 
size  n  =  100,  500  and  1000. 


ML  Estimator 

MM  Estimator 

n 

Pfa 

std|P/a 

Pfa 

Std  Pfa  (©>^)] 

100 

2.873E-4 

4.830E-4 

2.751E-4 

4.744E-4 

500 

1.286E-4 

0.961E-4 

1.280E-4 

0.985E-4 

1000 

1.141E-4 

0.582E-4 

1.138E-4 

0.596E-4 

The  false  alarm  performance  of  the  statistical  normalizer  assuming  Weibull  reverber¬ 
ation  in  the  presence  of  a  target-like  interference  is  examined  in  Figs.  7  and  8.  The 
PDF  of  log10  (Pfa  (©,  hj)  of  eq.  (40)  using  the  ML  parameter  estimates  is  shown  for 
the  statistical  normalizer  with  and  without  target  pruning.  Figure  7  represents  the 
case  of  n  =  100  auxiliary  data  samples  and  Fig.  8  when  n  =  500  samples  are  used. 
Both  have  5  auxiliary  data  samples  corrupted  by  a  non-fluctuating  (see  eq.  (16)) 
interference  with  the  indicated  interference-to-reverberation  ratio  (IRR)  and  the  re¬ 
maining  samples  unit  power  Rayleigh  distributed.  The  target  samples  were  inserted 
adjacent  to  each  other  with  the  position  chosen  randomly  for  each  of  the  10,000 
trials.  The  target  pruning  algorithm  assumed  a  maximum  of  mmax  =  10  target 
samples  when  n  =  100  and  7nmax  =  20  when  n  —  500  and  used  thresholds  generated 
according  to  eq.  (34)  with  a  =  0.1.  The  PDF  of  the  conditional  Pfa  when  there  is 
no  interference  present  is  either  the  upper  (n  =  100)  or  middle  (n  =  500)  graph  of 
Fig.  6.  As  the  IRR  increases,  the  PDF  of  the  conditional  Pfa  and  its  mean  for  the 
normalizer  without  target  pruning  shifts  far  to  the  left  of  the  desired  Pfa  =  10  4, 
indicating  much  higher  effective  thresholds  and  a  reduced  detection  performance. 
When  target  pruning  is  employed,  the  Pfa  does  not  change  substantially,  though 
the  conditional  Pfa  indicates  that  a  slight  reduction  in  detection  performance  will 
still  be  observed  as  there  is  more  mass  in  the  PDF  at  lower  values. 
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Figure  7  Probability  density  function  of  log10  of  instantaneous  Pfa  for  statistical 
normalizer  with  ML  parameter  estimator  with  and  without  NPIP  for  varying  IRR. 
Auxiliary  data  consisted  of  100  samples  of  Rayleigh  data ,  5  of  which  were  corrupted 
by  non-fluctuating  interference .  The  actual  Pja  for  each  detector  is  indicated  by  the 
*. 


Figure  8 


Same  as  the  above  figure  except  with  500  samples  of  auxiliary  data . 
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4.4  Statistical  normalizer:  Rayleigh  mixture  reverberation 

If  the  reverberation  amplitude  data  sample  X  is  distributed  according  to  a  Rayleigh 
mixture  with  parameters  0  =  [ni  •  •  •  ftp  \\  •  •  •  A ,p]  and  g  is  as  in  eq.  (14)  with 
the  parameter  estimates  0  =  [xi  •  ♦  •  7tp  Ai  •  •  •  Apj ,  then  the  instantaneous  Pfa 
is 

Pfa(&,h)  =^7riexpj--^  [g-1  (^;©)]2}  (41) 

where  g~l  (h;  ©)  is  the  functional  inverse  of  g,  which  does  not  have  an  analytical  so¬ 
lution  when  there  are  two  or  more  components  in  the  mixture  with  different  powers. 
If  the  maximum  likelihood  estimates  for  the  mixture  proportions  and  powers  are 
obtained  via  the  EM  algorithm  described  in  Section  2.2.2,  then  the  false  alarm  per¬ 
formance  is  scale  invariant;  that  is,  if  the  power  of  the  reverberation  data  changes, 
the  function  g  (x;  0^ ,  the  instantaneous  Pja  and  the  P/a  do  not  change1 . 

In  order  to  examine  the  performance  of  the  statistical  normalizer  in  realistic  rever¬ 
beration,  the  Rayleigh  mixture  parameters  are  estimated  from  the  SCARAB  data 
used  to  create  the  Pfa  estimates  of  Fig.  1  and  data  from  one  of  SACLANT  Undersea 
Research  Centre’s  SWAC  sea-trials.  In  both  cases,  the  data  were  well  fit  by  a  three 
component  Rayleigh  mixture  with  the  proportions  and  powers  as  noted  in  Table  2 
for  the  SWAC  data  and  Table  3  for  the  SCARAB  data.  Note  that  these  SWAC  data 
may  also  be  adequately  fit  by  a  two  component  Rayleigh  mixture.  An  indication  of 
the  severity  of  non-Rayleighness  of  these  distributions  is  seen  in  the  (/?2,/?i)  plots 
found  in  Fig.  9.  This  type  of  plot  shows  the  departure  of  the  third  and  fourth  order 
moments  (skewness  and  kurtosis)  from  that  of  the  Rayleigh  distribution  (see  [6]  for 
further  information).  Though  both  the  reverberation  models  exhibit  non-Rayleigh 
character,  the  SCARAB  data  model  is  extremely  non-Rayleigh  as  evidenced  by  the 
large  values  of  both  skewness  and  kurtosis. 

The  PDF  of  log10  of  the  instantaneous  Pfa  for  the  CA-CFAR  and  statistical  normal¬ 
izes  is  shown  in  Fig.  10  for  both  reverberation  models  when  the  desired  Pfa  is  10  . 

Clearly,  use  of  the  CA-CFAR  produces  a  substantially  larger  Pfa.  The  statistical 
normalizer  brings  the  Pfa  closer  to  the  design  value  at  the  expense  of  an  increase 
in  the  variance  of  the  logjQ  of  the  instantaneous  Pfa-  As  expected,  the  false  alarm 
performance  of  the  statistical  normalizer  becomes  closer  to  the  desired  performance 
as  more  data  is  used  to  estimate  the  Rayleigh  mixture  parameters.  That  the  sta¬ 
tistical  normalizer  is  not  invariant  to  changes  in  the  mixture  parameters  is  seen  in 
the  difference  in  performance  between  the  SWAC  and  SCARAB  data  models.  The 

1This  may  be  proven  by  first  showing  that  the  EM  algorithm  estimates  of  the  mixture  proportions 
{ifi}  are  not  affected  at  all  by  a  change  in  scale  and  then  noting  the  dependence  of  g  (x:  @)  on  the 
ratio  ^  which  becomes  distributionally  free  from  any  scale  applied  to  the  data. 
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less  extreme  reverberation  from  the  SWAC  data  model  is  easily  normalized  by  using 
500  samples  of  auxiliary  data  while  the  SCARAB  data  requires  more  than  twice  as 
much  to  provide  equivalent  control  over  the  false  alarm  performance. 

Figures  11  and  12  illustrate  how  the  false  alarm  performance  degrades  when  an  in¬ 
terference  5  samples  long  is  introduced  into  the  auxiliary  data  if  the  NPIP  algorithm 
is  not  first  used  to  identify  spatially  compact  interference  and  then  remove  it.  When 
the  NPIP  algorithm  is  not  used,  the  probability  of  false  alarm  is  substantially  lower, 
indicating  target  masking  and  a  loss  in  detection  performance.  Use  of  the  NPIP 
algorithm  results  in  minimal  change  in  performance. 

The  results  shown  in  Fig,  11  used  100  samples  of  auxiliary  data  while  the  results  of 
Fig.  12  used  500  samples.  The  number  of  order  statistics  used  in  the  NPIP  algorithm 
(f^rnax)  must  be  adjusted  according  to  the  amount  of  auxiliary  data;  in  particular, 
mmax  must  be  increased  when  more  auxiliary  data  are  used.  This  may  be  explained 
by  the  occurrence  of  a  greater  number  of  reverberation  data  that  exceed  the  inter¬ 
ference  data  values,  thus  requiring  consideration  of  a  larger  number  of  the  higher 
order  statistics  to  guarantee  capturing  interference.  In  the  simulations  performed 
for  this  report,  ramax  =  10  when  n  =  100  samples  of  auxiliary  data  were  used  and 
rnmax  =  20  when  n  =  500  samples  were  used. 


Table  2  Mixture  proportions  and  powers  estimated  from  SWAC  data . 


Proportion 

Power 

0.5958 

1.5529 

0.3288 

0.2262 

0.0754 

0.0053 

Table  3  Mixture  proportions  and  powers  estimated  from  SCARAB  data . 


Proportion 

Power 

0.4290 

0.9500 

0.4241 

0.6050 

0.1469 

2.2863 
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Figure  9  Skewness  and  kurtosis  plotted  on  a  (/%,  Pi)  plane  for  the  Rayleigh  mix¬ 
ture  models  estimated  from  the  SWAC  and  SCARAB  reverbemtion  data.  Note  the 
extremity  of  the  SCARAB  Rayleigh  mixture  model. 


SACLANTCEN  SR-303 


-  26  - 


SACLANTCEN  SR-303" 


SWAG  Data  Model 


C 

o 

o 

c 

3 

H— 

& 

(/} 

c 

CL) 

■o 

>N 


JQ 

CO 

n 

o 


n=1000 


Statistical 

normalizer 


/  \ 
t  \ 


CA-CFAR 


-4  -3.5  -3  -2.5 

log10  (Instantaneous  Pfa ) 


SCARAB  Data  Model 


Figure  10  Probability  density  function  of  log10  of  instantaneous  Pfa  for  statistical 
normalizer  and  CA-CFAR  detectors  with  100,  500  and  1000  samples  of  Rayleigh 
mixture  distributed  auxiliary  data.  The  actual  Pja  for  each  detector  is  indicated  by 
the  *. 
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Figure  11  Probability  density  function  of  log10  of  instantaneous  P/a  for  statistical 
normalizers  with  and  without  NPIP  for  varying  IRR.  Auxiliary  data  consisted  of  100 
samples  of  Rayleigh  mixture  data,  5  of  which  were  corrupted  by  a  non-fluctuating 
interference.  The  NPIP  assumed  a  maximum  of  10  corrupted  samples.  The  actual 
Pfa  for  each  detector  is  indicated  by  the  *. 
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Figure  12  Probability  density  function  of  log10  of  instantaneous  Pja  for  statistical 
normalizers  with  and  without  NPIP  for  varying  IRR.  Auxiliary  data  consisted  of  500 
samples  of  Rayleigh  mixture  data ,  5  of  which  were  corrupted  by  a  non-fluctuating 
interference.  The  NPIP  assumed  a  maximum  of  20  corrupted  samples.  The  actual 
Pfa  for  each  detector  is  indicated  by  the  *. 
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5 

Detection  performance 


Detection  performance  depends  directly  on  the  type  of  detector  used  and  the  sta¬ 
tistical  characterization  of  the  reverberation  and  target  echo.  To  understand  the 
relationships  between  detection  and  false  alarm  performance  in  the  presence  of  non- 
Rayleigh  reverberation,  first  assume  that  the  distribution  of  the  reverberation  is 
known  completely  and  that  a  simple  threshold  crossing  detector  is  utilized;  that  is,  a 
detection  is  declared  when  the  matched  filter  output  exceeds  a  certain  level.  The  the¬ 
oretical  receiver  operating  characteristic  (ROC)  curve  for  detecting  a  non-fluctuating 
or  Swerling  Type  0  target  in  Rayleigh  and  non-Rayleigh  reverberation  with  such  a 
detector  is  shown  in  Fig.  13.  The  non-Rayleigh  reverberation  has  the  same  power 
level  as  the  Rayleigh  reverberation  (thus  the  signal-to-reverberation  power  ratio  is 
the  same  for  both  curves)  and  is  in  the  form  of  a  Rayleigh  mixture  distribution  with 
the  parameters  as  found  in  Table  2.  If  the  desired  Pfa  is  10“4,  the  Pd  is  above  0.9 
for  Rayleigh  reverberation.  However,  if  the  data  are  non-Rayleigh  (with  the  distri¬ 
bution  previously  described)  and  a  perfect  mean  power  level  normalizer  is  used,  Pfa 
increases  by  over  an  order  of  magnitude  and  Pd  increases  very  slightly  (owing  to  the 
heavier  tails  of  the  reverberation),  as  illustrated  in  Fig.  13.  The  objective  of  the  sta¬ 
tistical  normalizer  is  to  account  for  the  distribution  of  the  reverberation,  essentially 
moving  along  the  dashed  curve  down  to  the  point  of  the  desired  Pja.  The  ensuing 
loss  in  Pd  represents  the  cost  of  controlling  the  false  alarms.  It  should  be  noted,  how¬ 
ever,  that  these  curves  assume  identical  signal-to-reverberation  power  ratios  (SRRs). 
This  may  not  be  the  case  in  practice;  for  example,  increasing  bandwidth  reduces  the 
range-bearing  resolution  cell  size,  potentially  creating  non-Rayleigh  reverberation, 
but  also  reducing  the  reverberation  power  level  (and  possibly  the  target  echo  level) 
resulting  in  a  different  SRR.  Fundamentally,  when  the  reverberation  distribution  is 
known,  the  detector  in  the  non-Rayleigh  reverberation  must  operate  somewhere  on 
the  dashed  curve. 

When  the  reverberation  distribution  is  not  known  exactly,  it  must  be  modeled  using 
an  appropriate  statistical  distribution  whose  parameters  are  then  estimated.  This 
leads  to  sources  of  errors:  mismatch  between  the  actual  reverberation  distribution 
and  the  assumed  one  and  losses  associated  with  the  estimation  of  the  parameters 
characterizing  the  reverberation  model.  Choosing  a  very  general  model,  such  as  a 
Rayleigh  mixture,  results  in  concentrating  the  losses  to  the  parameter  estimation. 
Choosing  a  specific  model,  such  as  the  Weibull  distribution,  may  lead  to  substantial 
mismatch  losses. 
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Figure  13  ROC  curve  for  non-fluctuating  target  in  Rayleigh  and  non-Rayleigh 
reverberation. 


The  conditional  Pfa  developed  in  the  previous  section  provides  a  technique  for  exam¬ 
ining  the  effects  of  parameter  estimation  and  the  losses  expected  from  distributional 
mismatch.  Given  enough  stationary  real  reverberation  data  is  should  also  be  feasible 
to  evaluate  the  conditional  Pfa  for  different  models  by  using  a  sample  CDF  in  place 
of  F(x\  0)  in  eq.  (39).  The  conditional  Pfa  can  also  provide  an  indication  of  detec¬ 
tion  performance  in  that  a  decrease  in  performance  is  expected  when  there  is  a  high 
probability  of  observing  values  of  the  instantaneous  Pfa  much  lower  than  the  desired 
Pfa.  This,  however,  is  not  concrete  enough  to  indicate  good  or  bad  performance. 

The  remainder  of  this  section  will  illustrate  through  simulation  that  it  is  feasible  to 
obtain  detection  performance  near  the  dashed  curve  of  Fig.  13  by  using  the  statis¬ 
tical  normalizer  of  Section  2  and  the  NPIP  algorithm  of  Section  3.  As  an  example, 
consider  a  non-fluctuating  target  occupying  5  adjacent  range  samples  and  auxiliary 
data  consisting  of  a  window  100  samples  long.  The  detector  in  this  case  simply 
compares  the  largest  value  of  the  target  bins  to  a  threshold.  The  probability  of 
detection  for  such  a  detector  is  shown  as  a  function  of  the  SRR  in  Fig.  14  when  the 
reverberation  is  Rayleigh  distributed  and  non-Rayleigh  using  the  Rayleigh  mixture 
of  Table  2.  Note  that  the  ordinate  has  been  plotted  using  the  standard  normal 
probability  scale  to  expand  the  regions  near  zero  and  one.  The  threshold  for  each 
situation  is  chosen  so  that  Pfa  =  0  .01  when  the  reverberation  is  Rayleigh,  as  is  re- 
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fleeted  in  the  curves  for  Rayleigh  reverberation  and  Rayleigh  mixture  reverberation 
where  the  statistical  normalizer  is  applied  with  perfect  knowledge  of  the  reverbera¬ 
tion  PDF.  It  can  be  seen  from  these  curves  that  the  desired  Pfa  is  achieved  by  the 
tendency  of  Pd  towards  this  value  as  SRR  decreases.  When  the  statistical  normalizer 
is  not  used  (i.e.,  assuming  the  reverberation  is  Rayleigh)  when  the  reverberation  is 
non-Rayleigh,  Pfa  exceeds  the  desired  value  along  with  an  increase  in  Pd.  The  cx’ 
(n  =  100)  and  ‘o’  (n  =  500)  marked  simulation  results  were  for  clean  auxiliary  data 
windows;  that  is,  the  data  were  all  distributed  according  to  the  Rayleigh  mixture  of 
Table  2.  These  curves  illustrate  that  increasing  the  amount  of  auxiliary  data  used  to 
estimate  the  reverberation  PDF  parameters  allows  achieving  performance  near  that 
when  the  PDF  is  known  exactly.  The  ‘+’  marked  simulation  results  placed  a  10  dB 
SRR,  5  range  sample,  non-fluctuating  interference  in  the  middle  of  the  auxiliary  data 
( n  =  100),  resulting  in  a  Rician  mixture  for  those  samples.  Here  it  is  seen  that  Pd  is 
slightly  below  the  performance  when  no  interference  exists,  indicating  that  the  NPIP 
algorithm  removed  most  of,  but  not  all  of,  the  interfering  signal  from  the  auxiliary 
data.  These  results  indicate  that  performance  near  the  expected  theoretical  values 
can  be  achieved  even  in  the  presence  of  interferences  by  using  the  NPIP  algorithm 
in  conjunction  with  a  statistical  normalizer,  with  improving  performance  as  more 
auxiliary  data  are  used. 
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Signal-to-reverberation  power  ratio  (dB) 


Figure  14  Probability  of  detection  as  a  function  of  signal-to-reverberation  power 
ratio  for  a  5  sample  target  with  a  detector  comparing  the  maximum  over  the  target 
cells  to  a  threshold.  The  thresholds  were  chosen  to  provide  a  false  alarm,  probability  of 
0.01.  The  statistical  normalizes  used  the  NPIP  algorithm  to  remove  any  interference 
present  in  the  auxiliary  data  prior  to  estimating  the  parameters  of  a  three  component 
Rayleigh  mixture.  The  V  (n  —  100 )  and  ‘o’  (n  —  500 J  marked  simulations  were  for 
clean  auxiliary  data;  that  is ,  auxiliary  data  with  no  interference  present.  The  ‘+’ 
(n  =  100 )  marked  simulation  had  auxiliary  data  contaminated  by  a  10  dB  SRR,  5 
range  sample,  non-fluctuating  interference . 
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6 

Application  to  real  data 


6.1  Experiment  and  sonar  configuration 


The  data  analyzed  in  this  report  were  collected  during  SACLANT  Undersea  Re¬ 
search  Centre’s  SCARAB  ‘97  sea-trial,  which  took  place  in  the  area  north  of  the 
island  of  Elba  off  the  coast  of  Italy.  Reverberation  data  were  collected  using  a 
high  power  towed  source  and  the  Centre’s  low  frequency  towed  line  array.  A  2  sec¬ 
ond  linear  frequency  modulated  (LFM)  waveform  sweeping  from  450  to  700  Hz  was 
transmitted  from  the  source,  which  was  towed  at  approximately  62  meters  depth 
and  was  comprised  of  three  transducers  spaced  approximately  1.5  meters  apart  and 
beamformed  to  point  40  degrees  down  from  horizontal.  The  hydrophone  data  were 
received  from  128  sensors  spaced  0.5  meters  apart  and  were  beamformed  using  Han¬ 
ning  shading  on  the  array  to  beams  spaced  so  that  their  beampatterns  overlapped 
3dB  down  from  the  main  response  axis  at  900  Hz,  resulting  in  54  beams  spanning 
from  forward  to  aft  endfire.  The  correspondence  between  beam  number  and  main 
response  axis  angle  is  depicted  in  Fig.  15  where  zero  degrees  indicates  forward  end- 
fire.  The  data  from  each  beam  were  matched  filtered,  basebanded  so  that  the  center 
of  the  waveform  band  (575  Hz)  shifted  to  zero  Hertz,  and  decimated  to  a  sampling 
frequency  of  250  Hz.  The  resulting  data  are  approximately  statistically  independent 
from  sample  to  sample  in  range  and  provide  a  range  sampling  interval  of  3  meters 
assuming  750  m/s  for  the  two-way  speed  of  sound. 


6.2  Single  beam  analysis 


One  of  the  beams  containing  non-Rayleigh  reverberation  is  chosen  for  more  detailed 
examination.  The  matched  filtered  beam  output  is  mean  power  level  normalized 
using  a  trimmed  mean  order  statistic  CFAR  processor  [2]  with  leading  and  lagging 
windows  each  containing  75  samples  with  a  buffer  of  2  samples  surrounding  the 
range  cell  being  processed.  The  order  statistic  CFAR  processor  trimmed  the  upper 
and  lower  10%  of  the  combined  leading  and  lagging  windows  of  auxiliary  data.  The 
configuration  of  the  normalizer  window  that  is  slid  along  the  matched  filtered  beam 
output  is  illustrated  in  Fig.  16. 
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Figure  15  Beamform,er  main  response  axis  angles  as  a  function  of  beam  number. 


Guard  Guard 

band  band 


Cells  being 
normalized 


Figure  16  Normalizer  configuration. 
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6.2.1  Distributional  models 

The  parameters  for  the  non-Rayleigh  models  described  in  [6]  are  estimated  and  used 
to  generate  frequency  curves  (CDFs  formed  by  using  estimated  rather  than  exact 
parameter  values),  from  which  the  probability  of  false  alarm  is  obtained,  as  shown 
in  Fig.  17.  Of  all  the  models,  the  three  component  Rayleigh  mixture  model  best  fits 
the  data  as  seen  in  Fig.  17  and  further  supported  by  the  p-values  of  the  Kolmogorov 
Smirnov  (KS)  test  statistics  between  the  models  and  the  data  found  in  Table  4  (a 
p- value  near  one  indicates  a  good  fit). 


Figure  17  Estimated  and  modeled  Pfa  for  5000  samples  of  data  from  the  middle 
region  of  beam  37.  The  Rayleigh  mixture  (3  component)  best  fits  the  data. 

Table  4  KS  p-values  for  various  non-Rayleigh  reverberation  models  with  data  from 
middle  region  of  beam  37. 


Model 

KS  p-value 

Rayleigh 

2.4e-14 

Log-normal 

3.3e-ll 

Edgeworth 

0.0014 

Weibull 

0.0082 

K 

0.3204 

Rayleigh  mixture 

0.9896 
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6.2.2  Statistical  normalization 

Statistical  normalizers  assuming  either  Weibull  or  Rayleigh  mixture  reverberation 
are  then  applied  to  the  mean  power  level  normalizer  output  under  the  assumption 
that  the  non-Rayleigh  character  of  the  reverberation  is  stationary  over  longer  ranges 
than  the  reverberation  power  level.  The  statistical  normalizers  are  implemented  us¬ 
ing  leading  and  lagging  windows  each  containing  250  samples  with  a  buffer  of  3 
samples  to  either  side  of  a  target  assumed  to  be  5  samples  wide  in  the  same  con¬ 
figuration  as  that  shown  in  Fig.  16.  To  decrease  the  amount  of  time  required  for 
processing,  the  same  nonlinearity  is  applied  to  all  5  of  the  assumed  target  sam¬ 
ples.  The  statistical  normalizer  assuming  Rayleigh  mixture  reverberation  always 
used  three  components  and  was  initialized  by  the  most  previous  mixture  parameter 
estimates  unless  one  of  the  mixture  proportions  or  powers  was  near  zero;  that  is,  if 
the  reverberation  were  well  fit  by  fewer  than  three  mixture  components.  The  data 
are  displayed  in  Fig.  18  along  with  the  threshold  that  produces  Pfa  =  10“4  in  unit 
power  Rayleigh  reverberation.  Clearly  the  statistical  normalizer  using  the  Rayleigh 
mixture  model  has  the  fewest  threshold  crossings  (false  alarms). 

The  ($2-  fix )  measures  of  the  skewness  and  kurtosis  [6]  are  estimated  from  overlapping 
blocks  of  1000  data  samples  for  the  mean  power  level  normalizer  and  the  statistical 
normalizers  assuming  Weibull  or  three  component  Rayleigh  mixture  reverberation. 
Figure  19  illustrates  how  non-Rayleigh  the  mean  power  level  normalizer  data  are,  and 
that  the  statistical  normalizer  using  a  three  component  Rayleigh  mixture  controls 
the  skewness  and  kurtosis  well. 

The  P/a  estimated  from  the  mean  power  level  normalizer  and  the  statistical  normal¬ 
izers  assuming  Weibull  or  three  component  Rayleigh  mixture  reverberation  is  shown 
in  Fig.  20.  Observe  that  the  statistical  normalizer  data  are  not  exactly  Rayleigh, 
though  their  false  alarm  performance  is  closer  to  the  desired  value  than  that  of  the 
mean  power  level  normalizer.  Were  the  reverberation  not  as  extremely  non-Rayleigh 
as  it  is,  as  discussed  in  Section  4.4,  the  false  alarm  performance  of  the  statistical 
normalizer  with  the  Rayleigh  mixture  model  would  be  closer  to  the  desired  Rayleigh 
reverberation  curve. 
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Figure  18  Normalized  data  for  beam  37.  The  dashed  line  represents  the  threshold 
for  which  the  Pfa  for  a  Rayleigh  random  variable  is  10-4. 


P2  =  kurtosis+3 

Figure  19  (/%,  Pi)  estimates  for  mean  power  level  normalizer  and  statistical  nor¬ 

malizes  assuming  Weibull  and  Rayleigh  mixture  reverbemtion.  The  estimates  were 
formed  using  blocks  1000  samples  long  from  the  data  in  beam  37  overlapping  50%. 
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Figure  20  Pja  estimated  from  beam  37  (approx.  10,000  samples )  for  mean  power 
level  normalizes  and  statistical  normalizers  assuming  Weibull  and  Rayleigh  mixture 
reverberation. 


6.3  Full  beam  analysis 


Images  of  all  the  beams  from  one  ping  of  data  are  shown  in  Fig.  21  in  a  rectangular 
range-beam  grid  for  the  mean  power  level  normalizer  and  the  statistical  normalizer 
using  a  three  component  Rayleigh  mixture  as  described  in  Section  6.2.2.  The  bright 
areas  visible  in  the  mean  power  level  normalizer  image  around  beams  10  and  35  are 
much  less  visible  in  the  statistical  normalizer  image.  Though  not  shown,  the  image 
for  the  statistical  normalizer  assuming  Weibull  reverberation  did  not  clean  up  the 
non-Rayleigh  reverberation  as  well  as  the  statistical  normalizer  using  the  Rayleigh 
mixture  model. 

To  provide  a  quantitative  measure  of  the  performance  of  the  statistical  normalizer, 
the  p- value  of  the  Kolmogorov- Smirnov  (KS)  test  statistic  between  the  observed  data 
(blocks  of  500  samples  with  90%  overlap)  and  the  Rayleigh  distribution  are  formed 
and  displayed  as  images  in  Fig.  22.  When  the  KS  statistic  p-valueis  near  zero  (black 
on  the  images),  the  Rayleigh  distribution  is  a  poor  fit  to  the  observed  data.  Values 
that  aren’t  small  indicate  that  the  data  may  be  well  fit  by  the  Rayleigh  distribution. 
For  example,  the  non-Rayleighness  of  the  direct  blast  (waveform  transmission)  is 
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Figure  21  Mean  power  level  normalizer  and  statistical  normalizer  data  using  a 
three  component  Rayleigh  mixture. 
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clearly  seen  in  the  forward  pointing  beams  (lower  beam  numbers)  for  all  of  the 
normalizes.  The  success  of  the  statistical  normalizer  in  converting  the  non-Rayleigh 
reverberation  observed  in  the  mean  power  level  normalizer  output  is  evident  in  the 
uniformity  of  the  KS  statistic  p-values  and  their  nearness  to  one.  Both  of  these 
indicate  that  the  statistical  normalizer  output  data  are  closer  than  the  mean  power 
level  normalizer  data  to  the  ideal:  unit  power  Rayleigh  reverberation  across  all  beam 
and  range  cells  where  ever  there  is  diffuse  reverberation. 

To  determine  how  close  the  normalizer  outputs  are  to  this  ideal,  a  similar  field  of 
Rayleigh  reverberation  was  generated.  The  KS  statistic  p- values  were  then  formed 
using  the  same  block  sizes  and  overlapping  as  for  the  normalizer  data.  Then,  as 
the  KS  statistic  p-values  were  themselves  random  variables,  the  CDF  of  the  p-values 
was  formed  from  the  data  beyond  4.5  km  (some  of  the  data  nearer  than  4.5  km  are 
invalid  owing  to  overloading  of  the  receive  array  hydrophones  by  the  transmission  of 
other  waveforms  in  the  wavetrain).  As  seen  in  Fig.  23,  the  CDF  of  the  KS  statistic 
p- values  for  the  mean  power  level  normalizer  indicates  a  much  higher  probability  of 
observing  non-Rayleigh  reverberation  (i.e.,  small  KS  statistic  p- values).  Application 
of  a  statistical  normalizer  assuming  Weibull  reverberation  reduces  this  probability 
somewhat,  but  not  as  much  as  when  the  Rayleigh  mixture  model  is  used.  The 
improvement  is  clearly  indicated  by  the  relative  nearness  of  the  Rayleigh  mixture 
curve  to  that  of  the  simulated  Rayleigh  data. 
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Mean  Power  Level  Normalizer 


Range  (km) 


Figure  22  KS  statistic  p-value  between  normalized  data  from  one  ping  after  mean 
power  level  normalization  and  statistical  normalization  assuming  Weibull  and  three 
component  Rayleigh  mixture  reverberation .  Values  near  one  more  closely  resemble 
the  Rayleigh  distribution . 
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Figure  23  Estimated  CDF  of  the  KS  statistic  p-values  of  Fig.  22  for  data  after 
the  first  4.5  kilometers. 
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7 

Conclusion 


Low  frequency  active  sonar  systems  operating  in  shallow  water  have  encountered 
non-Rayleigh  reverberation,  with  the  statistics  changing  in  range  and  bearing.  The 
lack  of  knowledge  about  the  reverberation  statistics  and  their  non-stationarity  hin¬ 
ders  the  application  of  automatic  detection  algorithms  to  range-bearing  images  ex¬ 
hibiting  such  non-Rayleigh  reverberation.  This  report  has  presented  a  method  for 
dealing  with  such  reverberation  by  converting  it  back  to  the  Rayleigh  distribution 
through  a  tail  de-emphasizing  non-linearity.  This  algorithm,  a  statistical  normal- 
izer,  was  evaluated  through  simulation  and  then  applied  to  real  data.  The  simu¬ 
lation  analysis  indicated  that  the  performance  was  linked  to  both  the  severity  of 
the  non-Rayleighness  of  the  reverberation  and  the  amount  of  data  used  to  estimate 
the  parameters  of  the  reverberation  model.  Application  to  real  data  indicated  that 
a  large  degree  of  the  non-stationarity  of  the  reverberation  could  be  removed,  pro¬ 
ducing  range-bearing  images  suitable  for  multi-beam  processing  in  the  presence  of 
non-Rayleigh  reverberation. 

In  applying  the  statistical  normalizer  to  real  data,  an  algorithm  was  developed  to 
prune  strong,  spatially  compact  interferences  from  auxiliary  data  that  are  subse¬ 
quently  used  to  estimate  reverberation  model  parameters.  The  algorithm  was  based 
on  the  indices  to  the  order  statistics,  and  is  therefore  non-parametric  (i.e.,  the  per¬ 
formance  is  independent  of  the  reverberation  distribution)  and  applicable  for  use  in 
reverberation  following  an  unknown  distribution. 

In  analyzing  the  false  alarm  performance  of  a  simple  threshold  crossing  detector 
operating  on  the  statistical  normalizer  output,  a  new  technique  of  examining  the 
effect  of  parameter  estimation  was  developed.  The  Pfa  conditioned  on  the  parameter 
estimate,  called  the  conditional  or  instantaneous  P/a,  was  used  to  determine  how 
variable  the  false  alarm  performance  can  be,  illuminating  the  effect  of  interference 
corrupting  auxiliary  data  and  the  effect  of  the  statistical  normalizer  on  non- Rayleigh 
reverberation. 

The  efficacy  of  the  statistical  normalizer  in  controlling  the  high  and  variable  false 
alarm  rates  exhibited  by  non-Rayleigh  reverberation  regions  observed  in  low-frequency 
active  sonar,  and  the  promising  simulation  results  indicating  acceptable  detection 
performance  lead  to  the  recommendation  that  this  algorithm  be  considered  for  in¬ 
clusion  in  the  Centre’s  real-time  low-frequency  active  sonar  system  for  use  and  the 
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eventual  evaluation  of  larger  data  sets.  At  the  same  time,  research  should  con¬ 
tinue  in  developing  reverberation  models  that  well  represent  real  data  with  a  more 
parsimonious  parameterization  than  the  Rayleigh  mixture  model. 
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Annex  A 

Target  models  in  Rayleigh  mixture  reverberation 


The  standard  fluctuating  and  non-fluctuating  target  models  may  be  represented  as 
special  cases  of  a  Swerling  Type  III  target  which  is  formed  as  an  additive  component 
(U)  to  the  complex  matched  filter  output  (Z), 

Y  =  Z  +  U  (Al) 


where  U  is  complex  Gaussian  distributed  with  mean  /x  and  variance  Ao  and  under 
the  Rayleigh  mixture  distributed  reverberation  assumption,  Z  is  distributed  as  a 
mixture  of  zero-mean  complex  Gaussian  random  variables  with  proportions  7 r*  and 
variances  A;  for  i  =  1,  . . . ,  p, 


v  1 

fz(z;Q)  =  yV*  — e 
»= 1  7rAi 


(A2) 


Here,  as  in  the  main  part  of  this  report,  0  =  [7Ti  •  •  ■  ttp  Ai  ■  •  •  Ap].  The  non¬ 
fluctuating  model  is  obtained  by  setting  Ao  =  0  and  the  fluctuating  model  by  setting 
(x  =  0. 


The  distribution  of  Y  may  be  found  by  first  conditioning  on  and  then  integrating 
over  C7, 


fY  (y;  0,  n,  A0) 


P 

5>< 

2=1 


7T  (Ao  +  Ai) 


(A3) 


which  is  seen  to  be  a  mixture  of  complex  Gaussian  random  variables  with  mean  (i 
and  variance  Ao  +  A*.  When  the  mixture  distribution  is  interpreted  as  one  of  the 
p  reverberation  models  occurring  with  its  associated  probability,  this  target-plus- 
reverberation  distribution  is  acceptable  on  an  intuitive  level  as  the  target  introduces 
the  same  effect  on  each  reverberation  component  when  it  occurs. 


By  transforming  the  real  and  imaginary  parts  of  Y  to  magnitude  and  angle  coor¬ 
dinates,  followed  by  integration  over  the  angle  coordinate,  the  distribution  of  the 
matched  filter  magnitude  output  ( X  —  |F|)  is  found, 


/(z;0,m,  Ao) 


2=1 

V 

x^ 


2x 

Ao  +  A  i 
2x 

Aq  +  A  i 


(x2  +  |^|2) 
.0  *0 +Ai 
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—  /l  (*>  X0  +  Ai,  fo)  (A-4) 

t=l 

where  &  =  J^2.  is  a  signal-to-reverberation  power  ratio  (SRR)  for  each  reverbera- 

Ao+Aj 

tion  component, 

fi  (*;  A,  5)  =  |e'{-T  J0  (2x^1  J  (A5) 

is  the  Rician  distribution  with  scale  parameter  A  and  SRR  parameter  6,  and  J0  (x)  is 
the  zero-order  modified  Bessel  function.  Thus,  the  magnitude  matched  filter  output 
is  a  mixture  of  Rician  random  variables  with  scales  Ao  +  A  j  and  SRRs  <%.  A  non¬ 
fluctuating  target  (Ao  =  0)  results  in  a  mixture  of  Rician  random  variables  with 

scales  \i  and  SRR  parameters 

/nf  (*;  M)  =  X) h  (x'>  Ai>  •  (A6) 

A  fluctuating  target  {p  =  0)  results  in  a  mixture  of  Rayleigh  random  variables  with 
power  Ao  +  Ai, 

p  9  rr  X 2 

fa  (®;  0,  Ao)  =  £><  *** .  (A7) 
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Annex  B 

Derivation  of  the  EM  parameter 
estimates  for  a  Rayleigh  mixture 


The  premise  of  the  EM  algorithm  lies  in  the  ability  to  describe  the  observed  data  as 
being  part  of  a  larger  set  of  data  with  a  simpler  likelihood  function  for  the  unknown 
parameters.  In  the  case  of  a  mixture  model,  there  is  a  very  natural  complete  set  of 
data.  If  a  random  variable  is  drawn  from  the  distribution  /  (2;  \)  with  probability 
7R  for  i  =  1,  . . . ,  p,  then  it  has  the  mixture  PDF  described  in  eq.  (12).  Observing  just 
the  reverberation  samples  {21,22,  •  •  •  ixn}  the  missing  data  are  simply  the  indices 
that  indicate  which  distribution  2 j  was  drawn  from;  that  is, 

Kj  =  i  (Bl) 

if  2 j  comes  from  mixture  component  i  for  j  =  1,  . . . ,  n.  The  PDF  of  the  complete 
data  is 

n 

f  (x,  K;  0)  =  JJ  7 TKjf  (xj]  \Kj)  (B2) 

3= 1 

where  x  =  [21  •  •  •  xn\T  and  K  =  [K\  •  •  •  Kn]T .  Note  that  Kj  can  take  on  integer 
values  from  1  to  p,  the  number  of  mixture  components.  Denote  as  the  set  of  all 
possible  values  of  K 

**  =  {1,  2,  (B3) 


The  PDF  of  the  missing  data  K  conditioned  on  the  observed  data  x  is 

/  (x,  K;  0) 


where 


/  (K  |x;  0)  = 


/  (x;  0) 

n"=i  VKjf  (xj;  \Kj) 

=  f[^„; 


/  (r;  ©)  =  n  E  ^  /  (®j;  Ai) 

j=1  i= 1 


(B4) 


(B5) 
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is  the  joint  PDF  of  the  observed  data  and 

WkjJ  =  H=1  n  f  ' 


(B6) 


The  expectation  step  of  the  EM-algorithm  requires  forming  the  expected  value  of 
the  log-likelihood  function  of  the  complete  data  where  the  expectation  is  taken  over 
the  missing  data  conditioned  on  the  observed  data, 


Q  (0,  0')  =  E©/  [log  /  (x,  K;  0)] 

=  £  log[/(x,K;0)]/(K|x;0') 


=  E  (t 

K e*K  U=1  L 

p  n 

=  EE4 

j'=l 


X  • 

log  7 TKj  +  log  2 Xj  -  log  A Kj  - 

xj 

log  7 Ti  +  log  2 Xj  -  log  Xi  —  — 


j=l 


(B7) 


where  W?-  is  formed  from  the  values  in  0'  and  the  simplification 


E  wt,  j  = 1 


(B8) 


has  been  used. 


The  maximization  step  of  the  EM-algorithm  requires  choosing  0  (i.e.,  the  Xj  and  A ») 
to  maximize  Q  (0, 0')  constraining  £Li  tt*  =  1  and  perhaps  requiring  each  individ¬ 
ual  proportion  7Tj  >  0.  Using  a  Lagrange  multiplier  to  enforce  the  first  constraint, 
it  is  easily  shown  that  maximizing  Q  (0,  0;)  with  respect  to  7 q  requires  choosing 


n 


1 


(B9) 


This  particular  choice  for  the  mixture  proportions  also  satisfies  the  latter  constraint 
mentioned  above  7 r*  >  0.  Maximizing  Q  (0,  0;)  with  respect  to  A,  results  in  choosing 


Ai  = 


£?=i  w'd  ■ 


(BIO) 
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Annex  C 

Mean  and  variance  of  Tm 


The  mean  and  variance  of  the  target  pruning  statistic  Tm  described  in  eq.  (29) 
must  be  determined  in  order  to  use  a  Gaussian  approximation  to  obtain  thresholds 
according  to  the  Type  I  error  probability.  Throughout  this  annex,  the  following 
summations  [29]  have  been  used  liberally, 


_  n(n  + 1) 

(Cl) 

^;2  n(rc+ 1)  (2n  + 1) 

h  "  « 

(C2) 

£,.3  n2  (n  +  l)2 
fc=i  4 

(C3) 

and 

^l4  n(n  +  1)  (2n  +  1)  (3n2  +  3n  —  1) 

h  ‘  “ 

(C4) 

It  is  instructive  to  first  consider  the  case  of  T<±  =  (i^i  —  A2)2  as  the  PMF  may  be 
obtained  analytically  and  used  to  verify  more  general  results.  As  K\  and  K2  can 
take  on  any  integer  values  from  1  to  n  as  long  as  they  are  not  equal,  there  are 
n  (n  —  1)  equally  likely  combinations.  Throughout  these  combinations,  there  are 
exactly  2  (n  —  j )  where  \K\  —  K2 1  =  j  for  j  =  1,  . . . ,  n  —  1.  Thus,  the  PMF  of  T2  is 

/2[f]  =  (4W  = 

1  0  otherwise 

(05) 

It  is  now  straightforward  to  determine  the  mean  T2, 

spy  - 

n  (n  +  1) 

6 

(06) 
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The  power  of  T2  may  be  similarly  obtained  through  direct  summation  over  the  PMF; 
however,  obtaining  it  from  the  definition  will  help  with  the  derivation  of  the  general 
case, 

E[T|]  -  E  [(Xx-l^)4] 

=  E  [K(\  -  4E  [KfK2]  +  6E  [k\K^  -  4E  [K^l]  +  E  [x24]  .  (C7) 

Due  to  the  symmetry  in  the  joint  PMF  of  the  indices  to  the  order  statistics,  the 
exact  index  or  indices  in  the  above  expectations  are  not  relevant.  For  example, 


ci  = 


E[Kt]=E[Kl]=E[Kf] 
n  .4 

= 


i—  1 


(n  +  1)  (2 n  +  1)  (3n2  +  3n  -  1) 
30 


(C8) 


and 


C2 


=  E  [KfK2\  =  E  [KxKl]  =  E  [KiK. 


1 


n  (n  —  1) 

1 


n  (  n 

EE*/-4 

i- 1  V”1 


i= 1  j= 1  i= 1 


4 


n(n  —  1) 

(n+1)  (I5n3  +  21n2  —  4) 


120 


(C9) 


In  the  latter  case,  note  that  using  summations  from  1  to  n  for  each  index  requires 
removing  the  terms  where  K\  =  K<2-  Similarly, 


C3 


1 

n  [n  —  1) 


1 

n  (n  —  1) 


E*4 


2=1 


(n  +  1)  (2 n  +  1) 
180  (n  —  1) 


(lOn3  —  3n2  — 


13n  +  6^  . 


(CIO) 


Substitution  into  eq.  (C7)  results  in 

E  [if]  =  2ci  -  8c2  +  6c3.  (Cll) 


SACLANTCEN  SR-303 


-54- 


SACLANTCEN  SR-303 


The  more  general  case  is  described  by  eq.  (29)  in  non-recursive  form  and  repeated 
here  for  convenience, 

m  i —  1 

(C12) 

i=2  j—\ 

The  mean  of  Tm  is  then  easily  found  from  the  mean  of  T2, 

E  [Tm]  = 


where  l  =  is  the  number  of  terms  in  the  double  summation.  The  power  of 

Tm  may  be  described  as 

m  i—  1  m  k— 1 

e  [ti]  =  e££Ee[(^-^)2(^-^)2]  (ci4) 

i=2  j=l  k=2  1=1 

=  Z  E  [(if!  -  K2f]  +  Zo  E  [(if !  -  K2)2  (if3  -  if4)2] 

+  [1(1-1)-  Zo]  E  [(ifi  -  if2)2  (ifi  -  if3)2]  ,  (C15) 

where 

‘-I?:;]-  <C16> 

The  four  summations  of  eq.  (C14)  result  in  l2  terms,  of  which  l  are  cases  where  the 
index  pairs  (i,j)  and  (fc,Z)  are  the  same  (i.e.,  i  =  k  and  j  =  l  or  i  =  l  and  j  =  fc), 
l  {l  —  l)  —  Zo  are  cases  where  the  pairs  share  only  one  common  index,  and  Zo  are  cases 
where  all  four  indices  are  different.  The  value  of  Zo  described  in  eq.  (C16)  is  obtained 
by  noting  that  for  each  set  of  four  indices  there  are  six  combinations  that  occur  in 
the  summations  of  eq.  (C14);  for  example,  the  combinations 

(if4-if3)2(if2-Kl)2,  (if4-if2)2(^3-^l)2,  (if4-ifl)2(if3-^2)2(C17) 

occur  as  do  their  reverse  multiplicative  order.  A  similar  argument  may  be  used  to 
show  that  the  number  of  combinations  with  one  common  index  is 

z(z-i)-Zo  =  eH  (Ci8) 


m,  i— 1 


i=2  j—1 
m  (m  —  1) 


m(m—  1) 


E[T2] 
n  (n  + 1) 


n  (n  +  1) 


6 


(C13) 
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The  third  expectation  results  in 

f  K\  -  2K\K3  +  K\K\  +  AKlK2K3  ' 

E  \{K1-K2)2{K1-Kzf\  =  E  -2KlK2  +  KlKl-2KlK1K2 

[  +K%Kl  -  2K%K1K3  _ 

—  ci  -  4 c2  +  3c3  +  0c5 

=  ci  —  4c2  +  3e3.  (C23) 

Substituting  eqs.  (C19),  (C20),  and  (C23)  into  eq.  (C15)  results  in 

e[t£]  =  l  (2d  -  8c2  +  6c3)  +  l0  (4c3  +  4c4  -  8c5) 

+  [l  (l  —  1)  —  Zo]  (ci  —  4c2  +  3c3)  .  (C24) 

The  variance  of  Tm  may  then  be  easily  obtained  from 

Var  [Tm]  =  E  fz£]  -  E  [Tm]2 .  (C25) 
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Annex  D 

Invariance  of  false  alarm  performance  in 
Weibull  reverberation  using  ML 
parameter  estimates 


If  the  maximum  likelihood  parameter  estimates  for  Weibull  distributed  reverbera¬ 
tion  are  used  in  the  statistical  normalizer  nonlinearity  of  eq.  (10),  the  false  alarm 
performance  is  invariant  over  the  whole  family  of  Weibull  distributions.  A  proof  of 
the  invariance  follows  in  this  annex. 


Either  by  differentiating  the  joint  PDF  of  Weibull  distributed  auxiliary  data  (Ai, 
. . . ,  Xn)  with  respect  to  a  and  (3  or  from  the  iteration  described  in  [6],  the  following 
relationships  may  be  derived  for  the  ML  parameter  estimates 

n  1  -1 


a  — 


1=1 


and 


L  _l±logXi  + 

K  "H  E”_,  A'? 


Letting  Zi  =  a  X?  these  relationships  may  be  described  as 


a  =  a? 


i  n  i 

n 


ir1 


(Dl) 


(D2) 


(D3) 


and 


ELi  Zi  log  Zi 
i 

En  rjP 

i= 1  * i 


m 


If  the  Xi  are  Weibull  with  parameters  a  and  /3,  then  the  Z{  are  exponentially  dis¬ 
tributed  with  unit  mean;  that  is,  the  Zi  are  distributionally  free  from  a  and  (3.  Thus, 
according  to  eq.  (D4),  the  ratio  |  is  also  distributionally  free  from  a  and  (3. 

Substituting  eq.  (D3)  into  the  nonlinearity  of  eq.  (10)  using  the  ML  parameter 
estimates  results  in 


g(X)  =  4&X& 
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z* 


Uzuzf\ 


(D5) 


where  Z  =  a  Xfi  is  exponentially  distributed  if  X  is  Weibull  distributed  with  pa¬ 
rameters  a  and  (3.  Thus,  g  (X)  is  distributionally  free  from  a  and  (3  resulting  in  a 
constant  false  alarm  probability  irrespective  of  the  actual  values  of  a  and  (3. 


Substituting  eq.  (D3)  into  the  instantaneous  Pfa  of  eq.  (40)  similarly  shows  a  lack 
of  distributional  dependence  on  a  and  /?, 


Pfa  (e,h) 


exp 


(D6) 
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